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Abstract 

We investigate the dynamics of a lattice soliton on a monatomic 
chain in the presence of damping and external forces. We consider 
Stokes and hydrodynamical damping. In the quasi-continuum limit 
the discrete system leads to a damped and forced Boussinesq equa- 
tion. By using a multiple-scale perturbation expansion up to second 
order in the framework of the quasi-continuum approach we derive 
a general expression for the first-order velocity correction which im- 
proves previous results. We compare the soliton position and shape 
predicted by the theory with simulations carried out on the level of the 
monatomic chain system as well as on the level of the quasi-continuum 
limit system. For this purpose we restrict ourselves to specific exam- 
ples, namely potentials with cubic and quartic anharmonicities as well 
as the truncated Morse potential, without taking into account exter- 
nal forces. For both types of damping we find a good agreement with 
the numerical simulations both for the soliton position and for the 
tail which appears at the rear of the soliton. Moreover we clarify why 
the quasi-continuum approximation is better in the hydrodynamical 
damping case than in the Stokes damping case. 

PACS: 63.10.+a Lattice dynamics: General theory and 05.45.Yv Soli- 
tons 



1 Introduction 

There is a long-standing interest in the dynamical and thermodynami- 
cal properties of anharmonic monatomic and diatomic chains ( see e.g. 
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|2|, |3], f|, ||, ^| ) . It was shown that these anharmonic chains can bear 
low-energy excitations which are solutions of a Boussinesq type equa- 
tion (in the long- wavelength approximation). For realistic interatomic 
potentials these soliton-like excitations are supersonic and correspond 
to a compression of the chain where the relation between amplitude 
(or width) and velocity depends on the form of the interatomic po- 
tential. They are very robust and propagate without energy loss, and 
their collisions are almost elastic even beyond the range of validity 
of the continuum approximation. Due to their robust character the 
soliton excitations are important in the coherent energy transfer and 
they have been used to explain energy transport in DNA j7j . There is 
also a growing evidence that nonlinear excitations participate in the 
heat conduction of anisotropic dielectric crystals §| || The 
non-diffusive heat flow was attributed to modified Korteweg-de-Vries 
solitons in 11 ] . The role of breathers for the thermal conductivity was 
studied in 0]. 

So far the main attention was paid to soliton dynamics in the ab- 
sence of dissipation. However, the dissipation influences significantly 
the solitons, changing their shape and velocity. An external driv- 
ing force is therefore necessary to sustain a soliton in steady state. 
The dynamics of slowly varying solitary wave solutions of the damped 
Korteweg-de-Vries equation was investigated by using the methods 
of inverse scattering theory |ll| , a multiple-scale perturbation expan- 



sion |jq| , and a Green's function formalism [17]. The soliton motion 
in Toda chains in the presence of dissipation and driving forces was 



studied in y, |lj. 



The objective of the present paper is to study properties of soli- 
tary waves in damped anharmonic lattices. We investigate two types 
of damping: Stokes friction (in Ref. [ 13 1 it is called outer friction) 
and hydro dynamical (or internal) friction []lg|j. We study the case 
of potentials with power-like anharmonicities as well as the case of 
a truncated Morse potential. We compare the results of numerical 
simulations, which are carried out for discrete anharmonic lattices as 
well as for the quasi-continuum Boussinesq system, with the results of 
a multiple-scale perturbation expansion obtained in the framework of 
the quasi-continuum approximation. 
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2 System and equations of motion 



We consider a chain of equally spaced particles of mass M (M = 1) 
with interatomic spacing a (a = 1) and displacement from equilibrium 
x n . The Lagrangian of our system is given by 

L = T-U- U ext . (1) 

Here 



T=\Y.^) (2) 

n 

is the kinetic energy ( x = =j^), 

U = Y J V(x n+1 -x n ) (3) 

n 

is the potential energy. We consider a potential between first neighbors 
of two types: the power-like potential 

V(r) = V harm + V anh , 

Vharm = \ r\ V anh = - 7* (p = 3, 4, ...) (4) 

2 p 
and the truncated Morse potential 

V(r) = -r 2 --r 3 + —r i (5) 
w 2 2 24 v ; 

which is the Taylor expansion of the Morse potential \ {e~ T - l) 2 . The 
last term in Eq. (|l|) represents the influence of external forces and is 
given by the expression 

Uext = ~ €n(t) {x n +l ~ X n ) (6) 
n 

To take into account damping effects in soliton dynamics we intro- 



duce the dissipation function ^[18|. In the case of the Stokes friction 



when the damping occurs due to interaction of the particles with a 
viscous environment (outer friction) the dissipation function depends 
on the velocity of the particles (Stokes law) and has the form 
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where v is the damping constant. In the case of internal friction 
(we will call this type of friction hydrodynamical), which is due to 
irreversible processes taking place within the system, the dissipation 
function depends on the time derivatives of the relative displacements 
and is given by the expression 

^h = ^^2(x n+1 (t)-x n (t)) 2 . (8) 

n 

The function (@) is the discrete version of the dissipation function 
which is usually used in macroscopic elasticity theory Jl8[] . 

It is necessary to point out that if we consider the soliton as a wave 
packet, its dynamics in the presence of Stokes damping shows that the 
long- wave components of the spectrum do not propagate (App. A). 
Thus the soliton decomposes after a time in the order of \/u. This 
feature does not show up with the hydrodynamical damping, if v < 1 
(see App. A). 

The equations of motion for the lattice displacements x n (t) in the 
presence of damping have the form 

d d L dL + _ Q 

dt dx n dx n dx n 

Substituting Eqs (|lj) - (||) into Eq. @, the equations of motion for 
the relative displacement u n = x n+ \ — x n can be written as 

u n = v'K+i) - + y'K„i) + 

tn+l{t)+in-l{t)-2t n {t) + D n (10) 

where V'(u) is the derivative of V with respect to its argument u and 
the damping term D n is determined by 

j—uiin for Stokes damping, 

, . . „ . « for hydrodynamical 

v (u n+ i + u n _i -2u n ) 
damping 

In order to obtain a analytical solution of the nonlinear system 
of equations (|l0|) we apply the quasi-continuum approximation pro- 
posed in |L9| (see also |2lJ] ) . Regarding n as a continuous variable 
(n — > x,u n (t) — > u(x,t)), Eq. (|To|) we obtain a damped and forced 
Boussinesq (Bq) equation (see App. B): 

d 2 t u - d 2 x u - d 2 t d 2 x u - dl (/(«)) = v m & x n d t u + d 2 Mx, t) 

(ii) 



4 



where d x and dt are the derivatives with respect to x and t, respec- 
tively; 



/(«) 



dV(u) 
du 



(12) 



is a nonlinear force and the right-hand-side of Eq. (11) represents 
the damping in the system and the action of an external force. The 
case m = corresponds to the Stokes damping while the case m = 2 
corresponds to the hydrodynamical damping: 



if m 
if m 



0, 
2. 



(13) 



3 Multiple scale expansion 

We are interested in how the dynamics and the behavior of the soliton 
is affected by the two types of damping (m = 0,2). So we consider 
both the position of the soliton center of mass as a function of time, 
X(t), and its shape for t > 0. We make a travelling wave ansatz 
u(x,t) = u(x — X(t)) and use a multiple-scale perturbation expansion, 
developed in detail in App. C, for a perturbed Bq equation 

d 2 u - d 2 u - d 2 d 2 x u - d 2 x (/(«)) = eF(x,t) (14) 

where eF(x,t) is the perturbation term with 

F{x,t) = v m d x n dtu + d?Z(x,t). (15) 

We seek an asymptotic solution of the form 

u = uq + e u\ + e 2 U2 H (16) 

with 

c = c + eci + e 2 c 2 H (17) 

where c = — <9t is the velocity of the soliton. e is a factor intro- 
duced for convenience in the analytical calculations. The case e = 
(u = uq) reduces Eq. (|14|) to the unperturbed Bq equation 

d 2 u - d 2 x u - d 2 d 2 x u - dl (f(u )) = , (18) 
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which is the well-known improved Boussinesq (IBq) equation 

] . When e = 1 we recover the damped and forced Bq equation ([ll]) . 
In order to interpret the multiple-scale perturbation results we must 
set e = 1 and assume that the terms on the r.h.s of (|ll|) are small 
enough. So, we must restrict ourselves to small values of the damping 
constant v. 

In what follows we restrict our study to the damped Bq equation 

d 2 u _ Q2 U _ Q2Q2 U _ Q2 {f{u)) = ^ d rn dtU (lg) 

The study of the effect of external forces, particularly stochastic forces, 
exceeds the frame of this paper and will be published later. 

From the multiple-scale perturbation analysis we obtain that there 
are two compatibility conditions: One of them follows from the order 
e 1 of perturbation, Eq. (p^); And the other one from the order e 2 , Eq. 
Both are valid for arbitrary potential V(u). 

Inserting the potential (Q) or (||), together with the corresponding 
one soliton solution, into the compatibility conditions we get a set of 
two ordinary differential equations of motion (ODEs). These ODEs 
govern the time evolution of the order e° and e 1 of velocity perturba- 
tion, namely cq and c\, respectively. 



3.1 The power-like anharmonic potential 



The expression (|63| ) is the one-soliton solution of the IBq equation 
flU) with the power- like anharmonic potential (Q). Substituting this 
solution in Eq. (p7|) and Eq. (|77f) yields 



(p-2)vc ( C Q 2 -1) 
6— 3p+2p cq 2 

c = { (20) 



(p+2)c (6-3p+2pc 2 ) 



m 



(p - 2) (3 (p - 2) + (18 - 7p) c 2 + 2pc 4 ) 

Cl = —V a Cl 

(6-3p + 2pc 2 ) 2 

+ v ^ ^ T(6 _ 3p + 2pC0 2 ) 4 r( _£_ )r( _jL_ ) 2 
x (6{p-2) 3 -2(p-2) 2 (-21 + 16p) c 2 
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+ p (104 - 122 p + 35 p 2 ) c 4 + p (16 + Up - 13p 2 ) c 6 

+ 2p 3 c 8 ) for m = (21) 



and 



^ Cl (p-2) 2 (c 2 -l) 2 
1_ (p + 2) 2 C 2 (6-3p + 2pc 2 ) 21 ^ P 4J 



+ V 



3 (-12 - 4p + 4 + 2p(p + 2) eg) 

2 2( P -2) 2 ^(cg-i) 3/2 r(|Ei) 2 r(^; 



(2 +p) 2 c 3 (6 _ 3p + 2pC0 2 ) 4 r( _£_ )r( _jL_ ) , 

3 (p - 2) 4 - 3 (p - 2) 3 (-10 + lip) cl 
+ (p-2) 2 p(-8 + 43p)^+p(-32 + 84p- 17p 3 )c(j 
+ 2p 3 (6+p)c[U for m = 2 (22) 



where T(-) is the gamma function. 

The set of ODEs (pO-p^) together with 



X = c where c = cq + c\ (23) 

constitute the complete set of ODEs which determine X(t) as a func- 
tion of time. 



3.2 The truncated Morse potential 

The one-soliton solution of the IBq with the truncated Morse potential 



is 



U °" l + J Bsinh 2 (2(x- Co t)) (24) 



where 



, _ ±6(c 2 - 1) B = JMj2lcl-12 



^3 + yj21($ - 12 ^3 + ^21c 2 - 12 

and 



4-1 
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(25) 



The upper sign means the soliton produces a rarefaction of the lattice, 
and the lower sign means compression. In this paper we consider 
only the compressional case. By inserting Eq. (24) in Eq. (pTj) and 



Eq. (77) we get a set of ODEs. They are rather cumbersome and 
therefore we make a further approximation by considering only the 
soliton dynamics close to the sound velocity. So we expand the ODEs 
in a Taylor series around the sound velocity where O(cq) terms are 
neglected. In the Stokes damping case the ODEs take the form 

68 19 2 49 

C0 + 45 I/C °- T5 UC °- 45" = ° 
2 v (-101 + 41 c ) 

Cl 15 (5 + 7 c ) Cl 

2v 2 y/2 (-39-1112c + 851cq 2 ) _ 
225VqT^T (5 + 7 c ) ~ ' 

and in the hydrodynamical damping case (m = 2) 

8 4 2 4 

co -— ^c + — vcq + — ^ = 
15 15 15 

, 32^(c -l) ^ 176 V2v 2 (cp-l)i 

Cl + , — r- Ci , r = U . (27) 

5 (5 + 7 c ) 75 (5 + 7c ) 

These approximations of the exact equations are good within a range 
of velocities 1 < cq ^1.1 



(26) 



Notice that either Eqs. (|26|) or Eqs. (|27|), together with Eq. (j2c 
constitute a complete set of equations of motion for X(t). 



4 Numerical simulations 

In order to verify these theoretical results for both the power-like and 
the Morse potentials we have performed molecular dynamics simula- 
tion for the discrete monatomic chain which is governed by Eq. flT0|). 
Moreover we have performed simulations on the level of the quasi- 
continuum limit, namely with the damped Bq equation (pi]). For the 
damped Bq system we have used finite-difference discretization in the 
space-domain |23| (see App. D for details). The time integration in 
both kinds of simulations was carried out by using the Heun method 
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p4|]. In order to start the simulations at t = we have used one- 
soliton solutions of the IBq equation. Since for some cases we need 
a long simulation time and the solitons are supersonic, we have used 
periodic boundary conditions to reduce the size of our system. In 
fact, we have used a chain with 1500 lattice points in the molecular 
dynamics simulation; and in the damped-Bq-simulation the length of 
the system has been L = 1000 with Ax = 0.25. We remark here that 
the solitons are bounded even when they develop a tail in the presence 
of perturbations, because the rear of this tail vanish. As the soliton 
including its tail is bounded in our finite system we are allowed to use 
periodic boundaries in our codes. The length of the tail grows with 
time, so we have considered not too long simulation times to avoid a 
possible overlapping between the rear of the tail and the front of the 
soliton. The other parameters have been At = 10~ 2 in the molecular 
dynamics simulation and At = 10 _1 in the damped-Bq-simulation. 

We have checked the accuracy of our codes by calculating the con- 
served quantity 

/oo 
u(x, t) dx 
-oo 

which is valid not only for the free soliton case {y m = 0) but also 
for the hydrodynamical damping case {m = 2). For the longest sim- 
ulation time the variation of this conserved quantity has been lower 
than 4 x 10 _9 % in the molecular dynamics simulation and lower than 
2 x 10 _13 % in the damped-Bq-simulation. Notice that this conserved 
quantity can be calculated in a numerical window as long as there is 
not overlap between its tail and the front of the soliton. 

The center of the soliton in both types of simulation has been 
found by finding the three points a?i_i, X{ and xi + \ where u{xi) is 
the absolute discrete maximum or minimum depending on whether 
the soliton is rarefactive or compressional. Afterwards a parabola has 
been fitted to the three coordinates, namely {xi,u(xi)}, and we have 
defined the vertex of this parabola as the soliton center of mass. 

We have chosen vq = —v = — 10 -3 and V2 = v = 10 -2 . The 
reason for choosing different values is that the Stokes damping has a 
stronger effect than the hydrodynamical damping for the same value v. 
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4.1 Soliton dynamics in the presence of hydro- 
dynamical damping (m = 2) 

In the case of the power-like anharmonic potential we have performed 
simulations in specific cases, namely for cubic as well as for quartic 
anharmonicity. We have also performed simulations in the case of the 
truncated Morse potential. 

In the cubic case (p = 3) in the presence of hydro dynamical damp- 
ing 

Eqs.d^H) reduce to 



C ° + 15(- C0 + 2c 3) 



-V 



rj 2 (-4-3T/ 2 + 2rj 4 + rf 
2 rf (22 + 40 ?7 2 + 18r/ 4 + rf 



Cl — V o ~ C\ 

15(l + r) 2 y 



225(1 + r, 2 f 

(28) 



where r\ has been defined in (p5[). In same way, for the quartic case 
(p = 4) Eqs. (|0H2|) reduce to 



. v (cq 2 - l) 2 

C0 + 3(-3c + 4 Co 3) 

rj 3 (-4 - 3?7 2 + 3r] A ) 

Cl — V Cl 

3(l + 3r/ 2 ) 2 
^ 2 7r 2 r? 5 (9 + 29r ? 2 + 39r ? 4 + 3r ? 6 ) ^ () 
~ V 36(l + 3?? 2 ) 4 ~ 



(29) 



In the case of the truncated Morse potential, which contains a combi- 
nation of cubic and quartic anharmonicities, we have already got the 
corresponding set of simplified Eqs. fl27|) . So depending on the type 
of anharmonicity we have solved either Eqs. (pa) or Eqs. (p^), or 



Eqs. (|27|) together with Eq. (23). The result from those numerical 
solutions is what we call theory's prediction. 

In Fig. |l] we show several examples of the soliton position as a 
function of time. These examples follow from the two kinds of simu- 
lations and from theory's prediction. In particular Fig. [l]a and Fig. 
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Figure 1: Soliton position in the sound velocity moving frame z vs. time 
in the hydrodynamical damping case with v = 10~ 2 . Solid line: molecular 
dynamics simulation, dashed line: theory's prediction (solution of Eqs. (|28|- 
|29D), dotted line: simulation for the Bq equation (0). Fig. (a) and (b) 
correspond to the cubic and quartic ^pharmonicities, respectively, and Fig. 
(c) corresponds to the truncated Morse potential. The uppercase letters A, 
B, C, D, and E correspond to different initial velocities, namely c(0) = 1.01, 
1.05, 1.1, 1.2 and 1.3, respectively 



[l]b correspond to the soliton dynamics of the cubic and quartic cases, 
respectively. And Fig. [l]c corresponds to the case of the truncated 
Morse potential. In each figure we show the soliton position for dif- 
ferent initial velocities. These cases are denoted by uppercase letters. 
The position is plotted in the sound velocity moving frame, defined by 
z(t) = X(t) — t. Notice that we consider lower initial velocities in the 
case of the truncated Morse potential (Fig. ||c), because in deriving 
Eqs. ( p7| ) further approximations have been made. In general for all 
the potentials we see that for low initial velocities (cases A, B and 
C) the lattice soliton position (solid-lines) is predicted rather well by 
theory's prediction (dashed lines) as well as by the position of the Bq- 
soliton (dotted lines). For higher initial velocities, namely c(0) = 1.2 
(cases D), the position of the lattice soliton agrees better with the 
position of the Bq soliton (dotted lines) than with the theory (dashed 
lines). And for even higher initial velocities, namely c(0) = 1.3 (cases 
E), quantitatively there is a clear difference between the three dy- 
namics. However qualitatively they are similar. The reason for this 
behavior is that the quasi-continuum approximation naturally predicts 
better the dynamics of broad solitons than that of narrow ones [^] ; 
and the solitons are narrower when the initial velocity is higher. No- 
tice that in the cubic and quartic cases theory(dashed lines) predicts 
better the Bq soliton position (dotted lines) than the lattice soliton 
position (solid-lines). This is due to the fact that both the theory 
and the Bq equation have been derived in the framework of the quasi- 
continuum limit. This feature does not show up in the case of the 
truncated Morse potential due to the further approximations that we 
have made in the theory of this case. 



4.2 Soliton dynamics in the presence of Stokes 
damping (m = 0) 

In this section we treat the same cases as in the previous section but 
with the soliton bearing systems in the presence of Stokes damping. 
As in the previous section Eqs. can be reduced depend- 



ing on whether the potential has cubic or quartic anharmonicity. In 
particular for the cubic anharmonicity (p = 3) 

, uc ( c 2 ~ !) n 
6 cq z — 3 
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+ + 2 -4 + 3?f + 8?f + 2r ? D 

Cl + Cl — 1/ o 1 = 0. 

3(l + r/ 2 ) 3 9r7(l + r7 2 ) 3 

(30) 

In the same way the reduced set of equations for the quartic case is 
v c (c 2 - 1) 

(r] + 3n 3 ) 2 (2- r? + 3r/ 4 ) 

Cl + 1/- 3 Cl — 

^(l + Sr/ 2 ) 4 
2 vr 2 (-1 - 3r? 2 - 5r? 4 + 19r? 6 +6r? 8 ) _ 
" 4r/(l + 3r ? 2 ) 4 

(31) 

And Eqs. (^7j) correspond to the case of the truncated Morse potential. 

The theory's prediction of the soliton position follows from the 
numerical solution of either Eqs.pO|) or Eqs. (31) or Eqs. ( |26| ) together 
with Eq. (P3|) . 

We show in Fig. ^ the same cases that we have treated in Fig. 
|]. We have also kept the same convention. In Fig. ||a we have 
not plotted the position of the Bq soliton because it agrees very well 
with the theoretical prediction (dashed line). In general we can ex- 
tend the comments made for the hydro dynamical damping case to the 
present Stokes damping case. We only want to remark some relevant 
differences. First, in the present case we have plotted the soliton dy- 
namics during the transient regime of the system, namely t < 1/v , 
in contrast to the hydro dynamical damping case where we have also 
considered times t >> \jv . It is because the overdamp character 
of the Stokes damping, in fact, the lattice soliton is destroyed by the 
damping, namely for times i^3/i/. On the other hand, neither the 
Bq simulation nor the analytical results predict in a correct way the 
lattice soliton dynamics for times t ^ 1/u . And second, the agreement 
between the position of the lattice soliton (solid-line) and either the 
position of the Bq soliton (dotted line) or theory's prediction (dashed 
line) is not as good as in the hydro dymamical case. For instance, if 
we compare the results for the quartic case with cq(0) = 1.2 (Figs, [ljb 
and |2|b: case D) we see that the agreement between simulations and 
theory is better for the hydrodynamical damping case than for the 
Stokes damping case. 
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Figure 2: Soliton position in the sound velocity moving frame z(t) vs. time 
in the Stokes damping case with v = 10~ 3 . Solid line: molecular dynam- 
ics simulation, dashed line: theoretical prediction (solution of Eqs. (|30l - [31~1) ), 
dotted line: simulation of the Bq equation (flTf). Fig. (a) and Fig. (b) cor- 
respond to the cubic and quartic anha^monicities respectively. And Fig. (c) 
corresponds to the truncated Morse potential. The uppercase letters A, B, 
C, D, and E correspond to different initial velocities, namely c(0) = 1.01, 
1.05, 1.1, 1.2 and 1.3, respectively 



4.3 Soliton profile. 



Up to now we have analyzed the dynamics of solitons but not their 
shape under the influence of damping. The soliton profile for t > 
in the presence of either hydrodynamical damping or Stokes damping 
can be obtained by a multiple-scale perturbation theory (App. C). In 
first order the soliton solution reads 

u = Uq + U\ (32) 

where no is the unperturbed solution (|^) and the function u\ follows 
from Eqs. ([79]) and (|S4|). As an example, the soliton solution in the 
case of a cubic potential with hydrodynamical damping reads 

1 

u = u + -M + w + v (33) 

with 

w = sech 2 (</>) (A x +A 2 (f) tanh(</>)) 

v = A 3 (ftsech 2 ((ft) + (At + (A 5 (ft 2 + A 6 )sech 2 ((/>)) tanh(</>) 
+ A 7 Tanh 3 {(j)) 

where (ft = r]9/2. M and the coefficients Aj, i = 1, .., 7 depend only on 
the time and are written down explicitly in App. E. In order to have 
a look at this theoretical behavior compared with the results from the 
simulations we show in Fig. ^ two specific examples which belong to 
the cubic case. Figs. ||a and ||b show snapshots of the soliton profile 
moving to the right side in the presence of hydrodynamical damping 
and Stokes damping, respectively. We see that in both examples the 
profiles are asymmetric and agree with each other rather well. The 
main feature which differs in both figures is that the amplitude of 
the tail which appears at the rear of the soliton in the presence of 
the hydrodynamical damping is positive, while in the other case it is 
negative. Notice that the theory's prediction of the Stokes damping 
case is not as good as in the hydrodynamical case where the dieviations 
are very small, in fact they are visible only in the center of the soliton. 

5 Discussion and conclusions 

In summary, in this work we have developed an analytical theory for 
the dynamics of lattice solitons on a monatomic chain under the in- 
fluence of damping. We have considered Stokes and hydrodynamical 
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Figure 3: Snapshots of the soliton profile for the system with cubic anhar- 
monicity. a: hydrodynamical damping case at t — 5000 with c(0) = 1.05, and 
v = 10~ 2 . b: Stokes damping case at t — 500 with c(0) = 1.1 and v = 10~ 3 . 
Solid circles: lattice soliton profile, dashed lines: Bq soliton profile, thin solid 
lines: theory's prediction. 
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damping. In the quasi-continuum approximation the dynamics of the 
driven and damped anharmonic lattice has been described by a driven 
and damped Boussinesq equation. Our analytical approach has been 
based on the multiple-scale perturbation theory. We have derived 
sets of equations of motion corresponding to zero- and first order per- 
turbations for the velocity. We have also calculated the first order 
perturbation for the soliton profile which develops a tail at the rear 
end. In order to check the validity of our results we have performed 
molecular dynamics simulation of the damped anharmonic lattice. We 
have also solved numerically the damped Bq equation. We have con- 
sidered lattices with cubic and quartic anharmonicities. The soliton 
position has been defined as the position of the maximum of the soli- 
ton. We have observed that our theory predicts in a correct way the 
dynamics of lattice solitons when they propagate in a medium with 
either hydrodynamical or Stokes damping. This good agreement also 
holds for the soliton profiles. However, in the Stokes damping case 
our analysis is only done for the transient time, namely t < 1/u where 
v is the damping constant, because the soliton decomposes for larger 
times. 

We have noticed that the quasi-continuum approximation describes 
in a better way the dynamics of the lattice solitons in the presence 
of hydrodynamical damping than in the case of the Stokes damping. 
This difference is due to the fact that in the hydrodynamical damping 
case the long- wave linear modes, which mostly contribute to the soli- 
ton dynamics, are underdamped (see Eq.(|43|)) while in the other case 
they are overdamped. 

In general, the agreement between our theory and molecular dy- 
namics simulation is mostly due to fact that our theory has been de- 
rived in the framework of the quasi-continuum limit. Moreover, this 
approach is better than earlier approximations made for the Korteweg- 
de Vries equation [16|, because higher soliton velocities can be consid- 
ered. 
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A Damping 



The Stokes damping does not permit the long wave components of a 
wave packet to propagate, while the hydro dynamical damping under 
a certain condition does not show this feature. This can be seen by 
means of a simple example: let us consider a harmonic monatomic 
chain with 2N lattice points whose equations of motion for the relative 
displacements u n (t) in the presence of Stokes damping have the form 

u n (t) = u n+ i(t) - 2u n (t) + u n -i(t) - uu n (t) (34) 

with n = 1,2, 3, 2N — 1. A travelling wave packet may be written 
as 

27V- 1 

u n {t) = u k e-^ n -^ (35) 

fc=0 

where (3 k is the wave number. 

By inserting (^) in (34) we get 2N equations of motion in /c-space: 

u? k - ivuik - 7fc = (36) 

where 

% = 2 (1 - <xx(p k )) (37) 
is a ^-dependent function which satisfies 

0<7fe<4 . (38) 



Solving (36) we obtain 



v I ( v x 2 



iuj k = --±id%- y-j . (39) 

Notice that for every finite value of v > there are small values of k 
for which the condition of oscillation 

7fc"(0 2 >O (40) 

is not satisfied. In other words, parts of the wave packet do not 
propagate. 
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However, this is not the case for the hydrodynamical damping. 
Here the equations of motion are 



u n {t) = u n+ i(t) - 2u n (t) + u n -i(t) 

+v (u n+ i{t) - 2u n (t) + u n -i(t)) 



(41) 



with n = 1, 2, 3, 2N — 1 By using the ansatz ([^) this set of equa- 
tions reads in A>space 



to k - iuj k uj k - 7 fc = 
with k = l,2,3,...,2iV - 1. 



\ 



Ik 1 



^ik 



Here the condition of oscillation 



1-1 g) ^ >0 



(42) 



(43) 



(44) 



is always satisfied, if v < 1 . 



B Quasi-continuum approximation 

The equations of motion of the lattice in the presence of dissipation 
and external forces are given by Eqs. (|To|) . In order to compact the 
calculations we consider here a dissipation term D n containing both 
the Stokes and the hydrodynamical damping. So in this case the 
equations of motion read 

U n = V'(u n+ l) ~ 2V'( Un ) + V'{u n -l) + 
£n+l(t) + U-l{t) ~2U(t) + VQUn + 

u 2 {u n +i - 2ii n + u n _i). (45) 

where u m , defined by Eq. (jl3|) , with m = 0, 2 are the damping 
constants of Stokes and hydrodynamical damping, respectively. This 
equation can be rewritten as 

u n = l{n) ( V'{u n ) + in{t) + u 2 u n ) + u u n (46) 



19 



where 

7(n) = 4 sinh 2 (^j = 4 sinh 2 {^d^j (47) 

is a differential operator where x = na and a is the lattice constant. At 
this point x is regarded as a continuous variable, so u n (t) — ► u(x, t) and 
£n(t) — * Taking into account that the function 4sinh 2 (<i/2)/d 2 

is smooth at d = 0, we can multiply both sides of (H) by the oper- 
ator a 2 9 2 /4sinh 2 (a5 x /2), and expanding this operator as well as the 
operator on the r.h.s. in a Taylor series we get 

d?u(x, t) = a 2 d 2 x V + a 2 d 2 x t(x, t) + a 2 \d 2 d?u{x, t) 
+U()dtu(x,t) — i>Qa 2 Xd 2 dtu(x,t) + U2a 2 d 2 dtu(x,t). 

(48) 

Setting a = 1, scaling x — ► y/\x, t — ► \/Ai, — * ^o/V^, ^2 — ► VX^2 
and using the definition (|l2|) we get 

9 2 n(x, t) - 9 2 n(x, t) - d 2 f(u(x, t)) - d 2 x d 2 u{x, t) = 

+<9 2 £(x, t) + u d t u(x, t) - i/ d 2 d t u(x, t) 

+u 2 dld t u(x,t). (49) 

One of the Stokes damping terms can be neglected because the field 
u(x, t) is slowly varying in space, therefore 

\uq d t u(x,t)\ >> \u d 2 d t u(x, t)\. (50) 

The estimate ( ]50D has been confirmed by the numerical solution of 
Eq.(|48|) with and without the term i>od 2 dtu(x , t) . In the rest of paper 
we consider separately either the Stokes damping case or the hydro- 
dynamical case, therefore Eq.(|49|) can be written as 

dtu - d 2 x u - d\d 2 x u - d 2 x (f(u)) = v m d^dm + d 2 x £{x, t) 

with m = 0, 2. 

C Multiple-scale perturbation expan- 
sion 

In this appendix we develop a multiple-scale perturbation approach 
to the generalized Boussinesq-Burgers equation 

dfu - d 2 x u - d 2 t d 2 x u - d 2 x (f(u)) = 
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e (v m d™d t u + d 2 x ti{x,t)) (51) 

where f(u) = dV d ^ — u is a nonlinear force and the right-hand-side of 
this equation represents the damping in the system and the action of 
an external force. We consider two types of damping 

{— v if m = 
v Mm = 2 

and e is a small parameter. In our derivation we will follow the proce- 
dure which was proposed in [jnj] for the perturbed Korteweg-de-Vries 
equation. 

By using the transformation to the moving frame of reference 
9 = x-X(T), T = et, 

T 

X{T) = i j c(T')dT' (52) 
o 

where X(T) is the center of mass position of the soliton and c(T) is 
its velocity which depend on the "slow" time variable T, Eq. (^T[) can 
be written in the form 

(c 2 dj -2ecd e d T - ecd e + e 2 o»|) (1 - dj) u 
-d$u-d 2 g (f(u)) = e(u m dZ l (ed T -cde) u + d 2 C(0,T)) 

(53) 

where ' = ^ and the notation £(x, t) = ((9,T) was used. We seek 
an asymptotic solution of the form 

u = u + e ui + e 2 u 2 H (54) 

with 

c = c + e ci + e 2 c 2 H (55) 



Inserting Eqs ( p4[ ) and ( |55[ ) into Eq. (^) and collecting powers of e 
we get 



a| ((eg - 1) u - cl d$ n - /(n )) = 0, (56) 
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e 1 : 



f 2 : 



d 2 e {{cl - 1) Ul - c 2 d% Ul - /'(n ) m) = d e F u (57) 

Fi = -(1 - dl) (2c c 1 dg - 2c d T - c ) n 

- z, m coCuo + <9 e C(#,T); (58) 

dl {(c 2 - 1) n 2 - cl 8$ u 2 - /'(n ) n 2 ) = dg F 2 
- (1 - <9|)«9t n + v m d™ (d T n - c <% ui - c\ dg n ) , 

(59) 

F 2 = -(1 - <9j){2 c ci <9 e ni + (c\ + 2c c 2 ) <% n - 
2 Co cV ni — 2ci<9t uo — cq U\ — c\ no} 

+\d e (f(uv)ut). (60) 

Integrating twice Eq. ( [56] ) under vanishing boundary conditions 
for no at infinity, we obtain the equation 

(eg - 1) n - cl dl u - f(u ) = 0. (61) 

In the case of the power-like anharmonic potential 

U anh = -uP, p>2 (62) 
P 



the solution of Eq. ( plj ) has the form 

1/(p-2) 



no 



(1(^-1)) ™ ' srf/<^>(^) (63 ) 

where the parameter t = , c ° characterizes the width of the excita- 
tion. 

Integrating Eq. (|57| ) under vanishing boundary conditions for dg u\ 
at 9 — > ±oo and for u\ — > as — ► oo we get 

oo 

(eg - 1) m - eg 9| ni - /'(n ) m = - J Fi(§) S. 

e 

(64) 
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The homogeneous part of Eq. ( p4[ ) has two linearly independent solu- 
tions 



vi = d g u , v 2 = v 1 



S 



{9§u ) 



(65) 



with the Wronskian v 2 d$ Vi — V\ d$ v 2 = 1. Taking into account that 
of these solutions only x>\ vanishes as 6 — > ±oo, the compatibility 
condition for a bounded solution of the inhomogeneous equation ( |64| ) 
is 



de Vl (6) Fx(e)de= / u Q F x de =o 



(66) 



Inserting Eq. (p8|) into Eq. (p6|) we obtain that the compatibility 
condition in the order e 1 has the form 



(67) 



d T ( c (uq + (d e u ) }) =v m c (u d™u ) + (u d 6 () 



where the notation 



(g)= / g{6)de 



(68) 



was introduced. In the same way, we obtain from Eq. ( p9[ ) that the 
compatibility condition in the order e 2 has the form 



Here 



2c («i (l - dfj u ) + ci (uq + (<% u ) 2 ) 

oo 

1 



(u (9) J d6d z T u Q (e)) + -{4 -l)M z = 
e 

v m (2co(«i <3^o) + ci (n Cn ) - (no cV S^uo) 



M= lim ui(9) 



(69) 
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and as seen from Eq. ( |57| ) this quantity is determined by the expres- 
sion 



M 



1 



cl-l 



(Fi(0)) 

2c dr{u ) + c (n ) - vq 5 m0 c (u Q ) 



+ C(oo,T)-C(-oo,r) 



(70) 



Note that in deriving the compatibility condition ( p9| ) the relation 

^(uod e (n?/"(«o))) = -( Ul d e { Ul f\u ))) (71) 
and Eqs (57)- (58) were used. 



Further simplification of the compatibility condition (69) may be 
achieved by using the relation 



(eg -!-/'(«„) -eg # 



duo 
del 



;i-$)«o 



(72) 



which can be obtained by differentiating Eq. (61) with respect to eg, 
and the relation 



1 



eg - 1 * "V dl 2 

which can be obtained from the equation 



it 2 dj - l) n + /(n ) = 



(73) 



(74) 



by differentiating it with respect to £ 2 . Thus by using Eqs ( pi| ) and 
© we get 

(»i (i - a|)u„) = - C1 — ((»§ + (a eui) ) 2 )) 

4c u <9c 

(75) 
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In the same way by using Eqs (JB1) and (f73|) we get 



(Ul 8g Uq) 







2c 



( U 2 - (deuof) + 



1 



1 



. V 2 ^cV ^C(e,T)--(C(-oo,T) + C(oo,r)) 



(76) 



Inserting Eqs ( |75| ) and (76) into Eq. (6S) we obtain that the compat- 
ibility condition in the order e 2 takes the form 



2 C ° 



d(u ) 



d 2 I 



dc 

i(C(-oo,T)-C(oc,r))^ 



V'm D m 



where the right-hand-side is determined by the expressions 



+ CiC 



2 c 



M (uq) - c 



d 
d c 



d(u ) \ 
dc J 
co 



+ ci <(cVo) 2 > 
d(u ) 



d c 



+ 



-a^o-^) fc(0,r)-i(C(-oo,r) 

co 0c / V 2 



+ C(oo,T)) 



+ 



D 2 

co 



<Tl M(« )-ci(^) 



2 c 



1 



d e u Q C(^T)--(C(-oo,r) + C(oo,T)) 



(77) 



(78) 



To find how the soliton profile changes in the presence of damping it 
is convenient to represent the function u\ in the form (see fi~6[| ) 

1 



u\ = -M + w + v 



(79) 
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where w — > as 9 — > ±00 while u — > T^-^ as # ^ ±oo. The functions 
^ and « satisfy the equations 

c 2 d 2 e w + (/'(u ) - eg + l) w = G w (6) (80) 
c 2 d 2 v + (/'(n ) - 4 + l) u = G°(9) (81) 



where 



G™(0) = 2c ci (1 - a|)n - ^Mf'(u Q ) (82) 



(7 

G v (e) = -J{(l- df) (2 c <9 T + c ) - u m c SP> u o (0) S 



(83) 

Using the functions ( |65| ) it is straightforward to see that the solutions 
of Eqs ( |8"0| ) and (^) are determined by the expressions 

e 

w = \ f (v 1 (e)v 2 (e)-v 2 (e)v 1 (e)) g w (9)S, 

9 

v = 1 j (y x {9) v 2 (9) - v 2 (9) Vl {9)) G v (9) d9. (84) 

Cq J 
u 

and ci to be obtained via Eq. (|77|). 



D Discretization of the Bq equation 

The Bq equations reads 

d?u(x, t) - 3%u(x, t) - dfd 2 x u{x, t) - d 2 x (f(u(x, t))) = 

F(x,t) (85) 

where F(x,t) are external forces and/or dissipation. By defining the 
variable v(x, t) = dtu(x, t) Eq. J85|) can be reduced to two partial 
differential equations of first order in time, namely 

dtv(x, t) = d 2 u(x, t) + d 2 dtv(x, t) + d 2 (f(u(x, t))) 
+F(x,t) 

dtu{x,t) = v(x,t). (86) 
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By using finite-difference discretization in the space-domain Eqs. fl86| ) 
take the form 

u i+ i(t) - 2ui(t) + Uj_i(t) 



+ 



+ 



Ax 2 

v i+ i(t) - 2vj(t) +Vj-i(t) 
Ax 2 

f(u i+1 (t)))-2f(u i (t)) + f(u i - 1 (t)) 
Ax 2 



where ' = 4, m(t) 



Ui(t) 

d_ 



+Fi(t), 

Vi{t) (87) 

u(xi,t), Vi(t) = v(xi,t), f(ui(t)) = u n (xi,t) and 
Fi{t) = F(xi,t) with n = 2,3. Xj = i Ax where Ax is the mesh size 
of the space variable and i = 1,2,- • -,N. The length of the system 
L = N Ax. In the numerical integration process we use periodic 
boundary conditions, namely uo(t) = ujv(^) and uw+i(t) = u\(t). 
The same boundaries are used for the variables Vi{t) and Fi{t). If we 
rewrite Eqs. (|87| ) so 

-« <+1 (t) + (Ax 2 + 2)vi(t) - Vi-i(t) = 
Ui+i(t) - 2ui(t) + Uj_i(i) + 
/(«i+l(t)) - 2f( Ui (t)) + /(«i_i(t)) + 
Ax 2 Fi(t), 

they can be regarded as a vectorial equations so 

Av = G, 



u 



where Ui and are elements of the vectors u and v, respectively. The 
elements Gi of the vector G are the r.h.s. of 
matrix 

(A -10 

-1 A -1 



(90) 
The 

and the square 
1 \ 











-1 A 



-1 








-1 











A -1 
-1 A J 



NxN 
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with A = Ax 2 +2. Notice that this tridiagonal matrix is cyclic because 
we use periodic boundary conditions [25|. From (|90|) we can derive 



v = A _1 G, 

u = v, (91) 

therefore at this stage we can use a classical integrator as for example 
the Heun algorithm in order to perform the numerical integration in 
time. 



E Coefficients 

The soliton solution in the case of cubic anharmonicity and hydrody- 
namical damping reads 

u = Uq + -M + w + v 

with 

w = sech 2 ((/>) (Ai + A 2 <j) tanh(^)) 

v = A 3 </>sech 2 (» + (A A + (A 5 cf 1 + A 6 )sech 2 (<fi)) tanh(^) 
+A 7 Tanh 3 (4>) 

where 



3 l /(3-c 2 ) v /^l 
Al = T 2cl -1 +3C1C ° 



3i/(3-cg)^/cg-l Cl 

A2 = 9 o — 

5 2 Cq — 1 Co 

A * = wS A 4 = -I(5cl-3)A 3 

A 5 = -^A 3 A 6 = -§(17 - lhcl)A 

A 7 = -l(5c 2 - 3)A 3 M = (5c 2 -3)A 3 
The velocities cq and c\ can be determined by Eqs. (| 



3 
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